Dyadic Filter Banks

using FileIO;
using LibSndFile;
using Plots;
using DSP;
using SampledSignals;
using FixedPointNumbers;
using CUDA
using Wavelets

tetris = load("tetris_piano_mono.wav");
fs = round(Int,tetris.samplerate);
tetris.data = tetris.data[1:(10*fs),:];
original = Float32.(tetris.data[:,1] .+ tetris.data[:,2]);
tetris

make_4d(x) = reshape(x,(size(x,1),1,1,1))
tensor(x) = CuArray(make_4d(x));
gpu_filt(h,x,p=(0,0),s=(0,0),d=(0,0)) = CUDA.CUDNN.cudnnConvolutionForward(
    h,
    x,
    padding = p,
    stride = s,
    dilation = d);
half(x) = size(x,1)÷2 + size(x,1)%2;
mag_resp(x) = 20.0 * log10.(abs.(freqz(PolynomialRatio(x,[1]))));

filters = [
    WT.makeqmfpair(wavelet(WT.haar)),
    WT.makereverseqmfpair(wavelet(WT.haar))
];
LA = tensor(Float32.(filters[1][1]));
HA = tensor(Float32.(filters[1][2]));
LS = tensor(Float32.(filters[2][1]));
HS = tensor(Float32.(filters[2][2]));

Analysis Filter Bank

x = [tensor(original)];
size.(x)

1-element Vector{NTuple{4, Int64}}:
 (441000, 1, 1, 1)

for depth ∈ 1:3
    y = cat(
        x[1:size(x,1)-1],
        [gpu_filt(HA, x[end], (1,0), 1, 1)[1:2:end-1,:,:,:]],
        [gpu_filt(LA, x[end], (1,0), 1, 1)[1:2:end-1,:,:,:]],
        dims=1)
    x = y;
end

size.(x)

4-element Vector{NTuple{4, Int64}}:
 (220500, 1, 1, 1)
 (110250, 1, 1, 1)
 (55125, 1, 1, 1)
 (55125, 1, 1, 1)

Synthesis Filter Bank

tmp = nothing;
y = nothing;
for octave ∈ 1:3
    y = cat(
        x[1:size(x,1)-2],
        [tensor(zeros(Float32,2*size(x[end],1)))],
        dims=1)
    tmp = tensor(zeros(Float32,2*size(x[end],1)));
    y[end][1:2:end] = x[end-1];
    tmp[1:2:end] = x[end];
    y[end] = gpu_filt(HS, y[end], (1,0), 1, 1)[2:end,:,:,:] .+ gpu_filt(LS, tmp, (1,0), 1, 1)[2:end,:,:,:];
    x = y;
end

R = Q0f15.(collect(x[1][:]));
recovered = copy(tetris);
recovered.data = hcat(R,R);
recovered